!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine weno_ly (fro, fmx, fmy, fmz, fbx, fby, fbz, fen, fpo,                                                                   &
   ro, mx, my, mz, bx, by, bz, en, po, nx, ny, nz, ch2)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny, nz

   real(8), intent(in) :: ch2
   real(8), dimension(nx, ny, nz), intent(in) :: ro, mx, my, mz, bx, by, bz, en, po
   real(8), dimension(nx, ny, nz), intent(inout) :: fro, fmx, fmy, fmz, fbx, fby, fbz, fen, fpo

   integer(4) :: jp1, jm1

   real(8), dimension(nx, ny, nz) :: ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, po_l
   real(8), dimension(nx, ny, nz) :: ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, po_r

   real(8) :: fro_l, fmx_l, fmy_l, fmz_l, fbx_l, fby_l, fbz_l, fen_l, fpo_l
   real(8) :: fro_r, fmx_r, fmy_r, fmz_r, fbx_r, fby_r, fbz_r, fen_r, fpo_r
   real(8) :: alpha_0, alpha_1, d_0, d_1, var_0, var_1
   real(8) :: vx, vy, vz, b2, v2, pr, c2, s2, ca2, cfx

!----------------------------------------------------------------------------------------------------------------------------------|
   d_0 = 2.0d0 / 3.0d0
   d_1 = 1.0d0 / 3.0d0
   do k = 1, nz
   do j = 2, ny - 1
      jp1 = j + 1
      jm1 = j - 1
      do i = 1, nx
         var_0 = 0.5d0 * (ro(i, j, k) + ro(i, jp1, k))
         var_1 = 0.5d0 * (3.0d0 * ro(i, j, k) - ro(i, jm1, k))
         alpha_0 = d_0 / ((eps + (ro(i, jp1, k) - ro(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (ro(i, j, k) - ro(i, jm1, k)) ** 2) ** 2)
         ro_l(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (mx(i, j, k) + mx(i, jp1, k))
         var_1 = 0.5d0 * (3.0d0 * mx(i, j, k) - mx(i, jm1, k))
         alpha_0 = d_0 / ((eps + (mx(i, jp1, k) - mx(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (mx(i, j, k) - mx(i, jm1, k)) ** 2) ** 2)
         mx_l(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (my(i, j, k) + my(i, jp1, k))
         var_1 = 0.5d0 * (3.0d0 * my(i, j, k) - my(i, jm1, k))
         alpha_0 = d_0 / ((eps + (my(i, jp1, k) - my(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (my(i, j, k) - my(i, jm1, k)) ** 2) ** 2)
         my_l(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (mz(i, j, k) + mz(i, jp1, k))
         var_1 = 0.5d0 * (3.0d0 * mz(i, j, k) - mz(i, jm1, k))
         alpha_0 = d_0 / ((eps + (mz(i, jp1, k) - mz(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (mz(i, j, k) - mz(i, jm1, k)) ** 2) ** 2)
         mz_l(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (bx(i, j, k) + bx(i, jp1, k))
         var_1 = 0.5d0 * (3.0d0 * bx(i, j, k) - bx(i, jm1, k))
         alpha_0 = d_0 / ((eps + (bx(i, jp1, k) - bx(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (bx(i, j, k) - bx(i, jm1, k)) ** 2) ** 2)
         bx_l(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (by(i, j, k) + by(i, jp1, k))
         var_1 = 0.5d0 * (3.0d0 * by(i, j, k) - by(i, jm1, k))
         alpha_0 = d_0 / ((eps + (by(i, jp1, k) - by(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (by(i, j, k) - by(i, jm1, k)) ** 2) ** 2)
         by_l(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (bz(i, j, k) + bz(i, jp1, k))
         var_1 = 0.5d0 * (3.0d0 * bz(i, j, k) - bz(i, jm1, k))
         alpha_0 = d_0 / ((eps + (bz(i, jp1, k) - bz(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (bz(i, j, k) - bz(i, jm1, k)) ** 2) ** 2)
         bz_l(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (en(i, j, k) + en(i, jp1, k))
         var_1 = 0.5d0 * (3.0d0 * en(i, j, k) - en(i, jm1, k))
         alpha_0 = d_0 / ((eps + (en(i, jp1, k) - en(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (en(i, j, k) - en(i, jm1, k)) ** 2) ** 2)
         en_l(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (po(i, j, k) + po(i, jp1, k))
         var_1 = 0.5d0 * (3.0d0 * po(i, j, k) - po(i, jm1, k))
         alpha_0 = d_0 / ((eps + (po(i, jp1, k) - po(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (po(i, j, k) - po(i, jm1, k)) ** 2) ** 2)
         po_l(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)
      enddo
   enddo
   enddo

   d_0 = 1.0d0 / 3.0d0
   d_1 = 2.0d0 / 3.0d0
   do k = 1, nz
   do j = 2, ny - 1
      jp1 = j + 1
      jm1 = j - 1
      do i = 1, nx
         var_0 = 0.5d0 * (3.0d0 * ro(i, j, k) - ro(i, jp1, k))
         var_1 = 0.5d0 * (ro(i, j, k) + ro(i, jm1, k))
         alpha_0 = d_0 / ((eps + (ro(i, jp1, k) - ro(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (ro(i, j, k) - ro(i, jm1, k)) ** 2) ** 2)
         ro_r(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (3.0d0 * mx(i, j, k) - mx(i, jp1, k))
         var_1 = 0.5d0 * (mx(i, j, k) + mx(i, jm1, k))
         alpha_0 = d_0 / ((eps + (mx(i, jp1, k) - mx(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (mx(i, j, k) - mx(i, jm1, k)) ** 2) ** 2)
         mx_r(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (3.0d0 * my(i, j, k) - my(i, jp1, k))
         var_1 = 0.5d0 * (my(i, j, k) + my(i, jm1, k))
         alpha_0 = d_0 / ((eps + (my(i, jp1, k) - my(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (my(i, j, k) - my(i, jm1, k)) ** 2) ** 2)
         my_r(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (3.0d0 * mz(i, j, k) - mz(i, jp1, k))
         var_1 = 0.5d0 * (mz(i, j, k) + mz(i, jm1, k))
         alpha_0 = d_0 / ((eps + (mz(i, jp1, k) - mz(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (mz(i, j, k) - mz(i, jm1, k)) ** 2) ** 2)
         mz_r(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (3.0d0 * bx(i, j, k) - bx(i, jp1, k))
         var_1 = 0.5d0 * (bx(i, j, k) + bx(i, jm1, k))
         alpha_0 = d_0 / ((eps + (bx(i, jp1, k) - bx(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (bx(i, j, k) - bx(i, jm1, k)) ** 2) ** 2)
         bx_r(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (3.0d0 * by(i, j, k) - by(i, jp1, k))
         var_1 = 0.5d0 * (by(i, j, k) + by(i, jm1, k))
         alpha_0 = d_0 / ((eps + (by(i, jp1, k) - by(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (by(i, j, k) - by(i, jm1, k)) ** 2) ** 2)
         by_r(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (3.0d0 * bz(i, j, k) - bz(i, jp1, k))
         var_1 = 0.5d0 * (bz(i, j, k) + bz(i, jm1, k))
         alpha_0 = d_0 / ((eps + (bz(i, jp1, k) - bz(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (bz(i, j, k) - bz(i, jm1, k)) ** 2) ** 2)
         bz_r(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (3.0d0 * en(i, j, k) - en(i, jp1, k))
         var_1 = 0.5d0 * (en(i, j, k) + en(i, jm1, k))
         alpha_0 = d_0 / ((eps + (en(i, jp1, k) - en(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (en(i, j, k) - en(i, jm1, k)) ** 2) ** 2)
         en_r(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)

         var_0 = 0.5d0 * (3.0d0 * po(i, j, k) - po(i, jp1, k))
         var_1 = 0.5d0 * (po(i, j, k) + po(i, jm1, k))
         alpha_0 = d_0 / ((eps + (po(i, jp1, k) - po(i, j, k)) ** 2) ** 2)
         alpha_1 = d_1 / ((eps + (po(i, j, k) - po(i, jm1, k)) ** 2) ** 2)
         po_r(i, j, k) = (alpha_0 * var_0 + alpha_1 * var_1) / (alpha_0 + alpha_1)
      enddo
   enddo
   enddo

   if (riemann_solver_flag .eq. 1) then
      call hlld_y (fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                         &
         ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, po_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, po_r, nx, ny, nz)
   elseif (riemann_solver_flag .eq. 2) then
      call hllc_y (fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                         &
         ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, po_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, po_r, nx, ny, nz)
   elseif (riemann_solver_flag .eq. 3) then
      call roe_y  (fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                         &
         ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, po_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, po_r, nx, ny, nz)
   endif

   do k = 1, nz
   do j = 1, ny - 1
      jp1 = j + 1
      do i = 1, nx
         vx = mx_l(i, j, k) / ro_l(i, j, k) 
         vy = my_l(i, j, k) / ro_l(i, j, k) 
         vz = mz_l(i, j, k) / ro_l(i, j, k) 
         b2 = bx_l(i, j, k) * bx_l(i, j, k) + by_l(i, j, k) * by_l(i, j, k) + bz_l(i, j, k) * bz_l(i, j, k)
         v2 = vx * vx + vy * vy + vz * vz
         pr = (en_l(i, j, k) - v2 * ro_l(i, j, k) / 2.0d0 - b2 / 2.0d0) * gmm1
         c2 = gm * pr
         s2 = c2 + b2
         ca2 = bx_l(i, j, k) * bx_l(i, j, k)
         cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_l(i, j, k) / 2.0d0)
         max_speed = cfx + abs (vx)

         if (riemann_solver_flag .eq. 0) then
            fro_l = mx_l(i, j, k)
            fmx_l = vx * mx_l(i, j, k) + pr + b2 / 2.0d0 - bx_l(i, j, k) * bx_l(i, j, k)
            fmy_l = vx * my_l(i, j, k) - bx_l(i, j, k) * by_l(i, j, k)
            fmz_l = vx * mz_l(i, j, k) - bx_l(i, j, k) * bz_l(i, j, k)
            fbx_l = po_l(i, j, k)
            fby_l = vx * by_l(i, j, k) - bx_l(i, j, k) * vy
            fbz_l = vx * bz_l(i, j, k) - bx_l(i, j, k) * vz
            fen_l = (en_l(i, j, k) + pr + b2 / 2.0d0) * vx -                                                                       &
               (bx_l(i, j, k) * vx + by_l(i, j, k) * vy + bz_l(i, j, k) * vz) * bx_l(i, j, k)
         endif
         fpo_l = ch2 * bx_l(i, j, k)

         vx = mx_r(i, jp1, k) / ro_r(i, jp1, k) 
         vy = my_r(i, jp1, k) / ro_r(i, jp1, k) 
         vz = mz_r(i, jp1, k) / ro_r(i, jp1, k) 
         b2 = bx_r(i, jp1, k) * bx_r(i, jp1, k) + by_r(i, jp1, k) * by_r(i, jp1, k) + bz_r(i, jp1, k) * bz_r(i, jp1, k)
         v2 = vx * vx + vy * vy + vz * vz
         pr = (en_r(i, jp1, k) - v2 * ro_r(i, jp1, k) / 2.0d0 - b2 / 2.0d0) * gmm1
         c2 = gm * pr
         s2 = c2 + b2
         ca2 = bx_r(i, jp1, k) * bx_r(i, jp1, k)
         cfx = sqrt((s2 + sqrt(s2 * s2 - 4.0d0 * c2 * ca2)) / ro_r(i, jp1, k) / 2.0d0)
         max_speed = max(cfx + abs (vx), max_speed)

         if (riemann_solver_flag .eq. 0) then
            fro_r = mx_r(i, jp1, k)
            fmx_r = vx * mx_r(i, jp1, k) + pr + b2 / 2.0d0 - bx_r(i, jp1, k) * bx_r(i, jp1, k)
            fmy_r = vx * my_r(i, jp1, k) - bx_r(i, jp1, k) * by_r(i, jp1, k)
            fmz_r = vx * mz_r(i, jp1, k) - bx_r(i, jp1, k) * bz_r(i, jp1, k)
            fbx_r = po_r(i, jp1, k)
            fby_r = vx * by_r(i, jp1, k) - bx_r(i, jp1, k) * vy
            fbz_r = vx * bz_r(i, jp1, k) - bx_r(i, jp1, k) * vz
            fen_r = (en_r(i, jp1, k) + pr + b2 / 2.0d0) * vx -                                                                     &
               (bx_r(i, jp1, k) * vx + by_r(i, jp1, k) * vy + bz_r(i, jp1, k) * vz) * bx_r(i, jp1, k)
         endif
         fpo_r = ch2 * bx_r(i, jp1, k)

         if (riemann_solver_flag .eq. 0) then
            fro(i, j, k) = 0.5d0 * (fro_l + fro_r - max_speed * (ro_r(i, jp1, k) - ro_l(i, j, k)))
            fmx(i, j, k) = 0.5d0 * (fmx_l + fmx_r - max_speed * (mx_r(i, jp1, k) - mx_l(i, j, k)))
            fmy(i, j, k) = 0.5d0 * (fmy_l + fmy_r - max_speed * (my_r(i, jp1, k) - my_l(i, j, k)))
            fmz(i, j, k) = 0.5d0 * (fmz_l + fmz_r - max_speed * (mz_r(i, jp1, k) - mz_l(i, j, k)))
            fbx(i, j, k) = 0.5d0 * (fbx_l + fbx_r - max_speed * (bx_r(i, jp1, k) - bx_l(i, j, k)))
            fby(i, j, k) = 0.5d0 * (fby_l + fby_r - max_speed * (by_r(i, jp1, k) - by_l(i, j, k)))
            fbz(i, j, k) = 0.5d0 * (fbz_l + fbz_r - max_speed * (bz_r(i, jp1, k) - bz_l(i, j, k)))
            fen(i, j, k) = 0.5d0 * (fen_l + fen_r - max_speed * (en_r(i, jp1, k) - en_l(i, j, k)))
         endif
         fpo(i, j, k) = 0.5d0 * (fpo_l + fpo_r - max_speed * (po_r(i, jp1, k) - po_l(i, j, k)))
      enddo
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine weno_ly
